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Abstract 


The upwind leapfrog or Linear Bicharacteristic Scheme (LBS) has previously been extended to treat 
lossy dielectric and lossy magnetic materials. This report extends the Linear Bicharacteristic Scheme for 
computational electromagnetics to the two-dimensional case, which includes treatment of lossy dielectric 
and magnetic materials and perfect electrical conductors. This is accomplished by implementing the LBS 
for homogeneous lossy dielectric and magnetic media and for perfect electrical conductors. Heterogeneous 
media are modeled by applying surface boundary conditions, and no special extrapolations or interpola- 
tions at dielectric material boundaries are required. The Perfectly Matched Layer (PML) outer boundary 
concept is also developed for this scheme. Results are presented for two-dimensional model problems on 
uniform grids, and the FDTD algorithm is chosen as a convenient reference algorithm for comparison. The 
results demonstrate that the explicit LBS is a dissipation-free, second-order accurate algorithm which uses 
an upwind computational stencil rather than a central difference stencil, and yet it has approximately one- 
third the phase velocity error. Computational requirements are also discussed. 


1 Introduction 

Numerical solutions of the Euler equations in Computational Fluid Dynamics (CFD) have illustrated the 
importance of treating a hyperbolic system of partial differential equations with the theory of characteristics 
and in an upwind manner (as opposed to symmetrically in space). These two features provide the motiva- 
tion to use the Linear Bicharacteristic Scheme ( LBS), also called the upwind leapfrog (UL) method, for the 
construction of many practical wave propagation algorithms. The upwind leapfrog (UL) method is based 
upon the Method of Characteristics, which is a widely used numerical solution concept in CFD [1]— [16]. 
In a hyperbolic system, the solutions (i.e. waves) propagate in preferred directions called characteristics. A 
characteristic can be defined as a propagation path along which a physical disturbance is propagated [17], 
The relevance to Maxwell’s equations is intuitively obvious because electromagnetic waves have preferred 
directions of propagation and finite propagation speeds. Characteristic-based methods have also been suc- 
cessfully implemented and demonstrated primarily for free space and perfect electrical conductor (PEC) 
electromagnetic problems [ 1 8]— [3 1 ] . 

This report extends the LBS to the two-dimensional case to model both homogeneous and heterogeneous 
lossy dielectric and magnetic materials and perfect electrical conductors (PECs). The LBS was originally 
developed to improve unsteady solutions in computational acoustics and aeroacoustics [32]-[38], It is a clas- 
sical leapfrog algorithm, but it uses a one-sided (or upwind) stencil for the spatial derivatives, which follows 
the wave characteristic more closely when compared with a classical leapfrog method. This approach pre- 
serves the time-reversibility of the leapfrog algorithm, which results in no dissipation, and it permits more 
flexibility by the ability to adopt a characteristic based method. Clustering the stencil around the character- 
istic enables high accuracy to be achieved with a low operation count in a fully discrete way [33]. The use 
of characteristic variables allows the LBS to treat the outer computational boundaries naturally using the 
exact compatibility equations. The LBS treats the outer boundary condition naturally without nonreflecting 
approximations. The interior point algorithm predicts the outgoing characteristic variables at the domain 
boundaries. For multidimensional applications, in principle, through knowledge of the wave propagation 
angle, the local coordinates can be rotated to align with the characteristics, at which the boundary condition 
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becomes almost exact. Therefore, no extraneous boundary condition is required. In the cases where this 
coordinate transformation is not implemented, the characteristic-based algorithm provides only an approxi- 
mation at the outer grid boundaries. However, the Perfectly Matched Layer (PML) outer boundary concept 
can be applied to this scheme, which is discussed later in this report. The LBS also offers a natural treatment 
of dielectric interfaces, without any extrapolation or interpolation of fields or material properties near mate- 
rial discontinuities. Exact boundary conditions on the tangential field components are directly enforced at 
material interfaces. The LBS offers a central storage approach with lower dispersion than the Yee algorithm 
[39], It has previously been applied to two and three-dimensional free-space electromagnetic propagation 
and scattering problems [34], [37], and it was recently extended to treat lossy dielectric and magnetic mate- 
rials for the one-dimensional case [40], 

The objective of this report is to present the extension of the LBS to the two-dimensional case, which 
includes lossy dielectric and magnetic materials. Results are presented for several two-dimensional model 
problems, and the FDTD algorithm is chosen as a convenient reference for comparison. The principles to 
extend this procedure to the three-dimensional case are straightforward. Sections 3 and 4 present the LBS 
implementation for the TM and TE polarizations, respectively. Section 5 outlines the dielectric material 
surface boundary condition and Section 6 discusses the outer radiation and PML boundary conditions. Sec- 
tion 8 reviews the Fourier analysis and computational requirements. Finally, Section 9 presents results for 
two-dimensional model problems and Section 10 provides concluding remarks. 

2 Abbreviation List 

The following table provides a list of abbreviations and acronyms used throughout this report. 

Abbreviation Description 


CFD 

Computational Fluid Dynamics 

FDTD 

Finite Difference Time Domain 

LBS 

Linear Bicharacteristic Scheme 

PEC 

Perfect Electrical Conductor 

PML 

Perfectly Matched Layer 

TE 

Transverse Electric 

TM 

Transverse Magnetic 

UL 

Upwind Leapfrog 

2D 

Two-dimensional 


3 TM Polarization 


Maxwell’s equations for linear, homogeneous and lossy media in the two-dimensional TM case (taking 
d/dz — 0) are 
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where a and a* are the electric and magnetic conductivities, respectively. Using the electric displacement 
D — eE and making the substitution c — 1 /^/ye gives 
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The procedure for the LBS is to transform the dependent variables D z , H x and H y to characteristic variables. 
The algorithm developed here is the simplest leapfrog scheme described by Iserles [41] combined with 
upwind bias, or simply, the Linear Bicharacteristic Scheme (LBS). To transform (4)-(6) into characteristic 
form, we multiply (5) and (6) by c and then add and subtract from (4) to give 
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Note that these equations are almost identical to the equations for the one-dimensional case [40], except for 
the addition of the cross-derivative magnetic field terms. The characteristic variables are defined as 


P = D z — - Hy 

c 

Q = D z + -Hy 

c 

R = D z + - H x 
c 

S = D z — - H x 
c 


(11) 

(12) 

(13) 

(14) 


to represent the ±x and ±y propagating solutions, respectively. Using these definitions, (7)— (1 0) can be 
rewritten as 
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It is convenient to define and store the following coefficients before time-stepping begins 
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Equations (15)— (18) can be rewritten more concisely as 
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To develop the discretized algorithm for a two-dimensional system, the stencils of Figures 1 and 2 are 
proposed for the LBS. We discretize time and space as t — nAt, x — iAx, y — jAy. To solve the wave 


(a) 


(b) 



Figure 1 : Two-dimensional upwind leapfrog computational stencils for right-going (a) and left-going (b) x 
propagating characteristics. 

propagation problem without introducing dissipation, it is necessary that the stencil have central symmetry 
so the scheme employed is reversible in time [33]. The stencil in Figure 1 a is used for Px propagating waves 
and the stencil in Figure lb is used for —x propagating waves. The upwind bias nature of these stencils is 
clearly evident. Figures 2a and 2b show the stencils for ±y propagating waves, respectively. References 
[32], [33], [36], [37], [38] clearly show that the LBS is second-order accurate. 

Note that the third and fourth terms in (2 1 )— (24 ) represent the electric and magnetic loss (or source) 
terms. A key element in developing an accurate LBS scheme is proper treatment of these source terms. 
The method used here indexes the self source term in (21) (i.e. P) at time level n + 1 and it indexes the 
coupled source term Q at time level n. This avoids a matrix solution at each grid point, and the formulation 
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Figure 2: Two-dimensional upwind leapfrog computational stencils for right -going (a) and left -going (b) y 
propagating characteristics. 

easily limits to the perfect conductor condition as a ->• oo. An identical application is made for equations 
(22)-(24). 


Using the stencils shown in Figures 1 and 2 and the source term indexing scheme described above, the 
resulting finite difference equations for (2 1 )-( 24 ) are 
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where i+" denotes the value for P at grid point (i, j) and time level n. Note that the differences are taken 
with respect to the cell center, i.e. the coordinate (i.j) is located at the center of the cell. Since we know 
that H x — c(R - S )/ 2 and H y — c(Q - P)j 2, these equations can be rearranged in the form 

(l+aAt)P+\} 2>j = P+^ + (l-2 Vx )(p+^-P+, 2 ^-bAtC+ +l ,^- 
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where — cAt/Ax and v y = cAt/Ay are the x and y Courant numbers. We now rewrite equations 
( 29 )-( 32 ) as 
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Equations ( 33 )-( 36 ) are the update equations for the 2 D TM LBS scheme at cell (*, j) which can contain 
lossy dielectric and magnetic materials. Note that as ct -> 00, then we have the PEC condition that 77 / 2 ? > 

^"-1/2,7’ 777/2’ and 777/2 - 0 as required. 

4 TE Polarization 

Maxwell’s equations for linear, homogeneous and lossy media in the two-dimensional TE case (taking 
d/dz — 0) are 
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Using the electric displacement D — eE and making the substitution c — 1 /y/JIe gives 
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The procedure for the LBS is to transform the dependent variables Z) x , and H z to characteristic variables. 
To transform (44)-(46) into characteristic form, we multiply (46) by c and then add and subtract from (44) 
and (45) to give 
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The characteristic variables are defined as 


+ -D x + c d fL + °—h z = 0 
ax /ic 


P = Py + ~H Z 

C 


Q = Py~~H z 
c 


R = P x --H z 
c 

S — P x 4 — H z 
c 


to represent the ±x and ± y right and left propagating solutions, respectively. Using these definitions, 
(47)— (50) can be rewritten as 
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Using the a and b coefficients defined in (19) and (20) we can rewrite equations (55)-(58) more simply as 
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To develop the discretized algorithm for a two-dimensional TE system, we use the same stencils as for the 
TM case, which are shown in Figures 1 and 2. We also employ the same indexing scheme for the self and 
coupled source terms in (59)-(62) and we also use a central difference approximation at the appropriate 
half-integer indexed cell to evaluate the cross derivative terms. 

To derive the finite difference equations for (59)-(62) we use the same stencils shown in Figures 1 and 
2. Since we also know that D x — (R + S) /2 and D y — (P + Q) /2, the TE finite difference equations are 
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(72) 
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Equations (67)-(70) are the update equations for the 2D TE LBS scheme at cell (?, j) which can contain 
lossy dielectric and magnetic materials. Note that as a -» oo, then we have the PEC condition that tf+ l 1 / 2 } -, 
Q”_jw 2 ? , and -S’”j l I 1 1 /2 = 0 as required. Note that the update equations are identical to the TM 

case, the differences being in the definition of the characteristic variables and in evaluation of the cross 
derivative terms. 

5 Heterogeneous Materials 

One of the difficulties with the conventional FDTD algorithm is the error in treatment of material discon- 
tinuities. Recent research efforts have attempted to reduce this error source by suitable averaging of material 
properties across the interface or by interpolation or extrapolation of the electromagnetic fields near these 
material boundaries [42], [43]. The advantage of the LBS is that the characteristic based nature of the al- 
gorithm leads to a very natural treatment of dieletric interfaces. Since the LBS works with characteristic 
variables, the slope of characteristic curves in each material will be different, and the physical boundary 
conditions permit an elegant and efficient implementation of a dielectric interface boundary condition. This 
numerical boundary condition implements the physics exactly, with no averaging, interpolation or extrapo- 
lation required. 

To implement the dielectric material interface boundary condition, consider a portion of a two-dimensional 
grid shown in Figure 3, which contains material discontinuities in both the x and y directions. We can see 
that the characteristic variables P and Q are co-located at the center of the cell edges along the y axis. Sim- 
ilarly, variables R and S are co-located at the center of the cell edges along the x axis. Thus, the LBS has 
a staggered storage scheme, similar to the conventional FDTD method. Spatial derivatives are taken with 
respect to the cell center, which is where the cell coordinates (i,j) are defined. 

The characteristic variables at each grid point (i,j) on the interface are split into two components each: 
P‘ 2 ,j and Q 2 ,? for interfaces perpendicular to the x axis and R t j , ll,^, and S-i- 2 for interfaces 
perpendicular to the y axis. The terms P\. j, Qi,j, Ri, 1 and S h \ exist just to the left and bottom of the ma- 
terial interface, respectively, as shown in Figure 3. The remaining terms Q 2 ,j, Ri , 2 and S ) : 2 exist just 
to the right and top of the material interface. Note that the i and j subscripts have been omitted from the 
dielectric boundary split field components in Figure 3 for clarity. For material 1, equation (33) is used to 
predict the value for T^T 1 at the boundary and for material 2, equation (34) is used to predict the value for 
' ■ Similarly, equation (35) is used to predict the value of /^’j hl and (36) predicts the value for S”^ 1 . 
The procedure for the TE polarization is identical. For example, in the TE case, the characteristic variable 
P uses field components D y and H z , which both are tangential to material interfaces that are perpendicular 
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Figure 3: Section of a two-dimensional computational grid for the LBS showing characteristic variables, 
dielectric interfaces and corresponding field components and characteristic variables used for the surface 
boundary condition. 

to the x axis. 

To complete the implementation, the Q 1 ^ 1 and P"J l terms must be updated. These terms are updated 
by enforcing the physical boundary conditions on the electromagnetic field at the material boundary. We can 
then solve for Q "+ 1 and P$+ l in terms of the “known” characteristic variables and To develop 
this procedure, the electromagnetic boundary conditions on the tangential field components are given by 

E zl .i - E z 2, (75) 


For the right-going wave, substituting (75) and (76) into (1 1) gives 
r?rc-fi n"+i -L ^ rr^+i 

P ,3 ~ + Cl 'V,; 

— — ( P n+1 + 0 n+1 'i + — ( } 
~ 2e 2 [ ^ ^ j + 2 Cl v 

Similarly, substituting (75) and (76) into (12) yields 


Q& 1 ) 


= D ” 


e ' 2 ( pn+l 

26 , \ ^ 


+ err 1 




Since P"^ 1 and are determined at boundary point (i,j) from the usual update equations (we treat them 
as “known” variables), it is necessary to express P, 1 ^ 1 and 1 in terms of these variables. Rearranging 
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(78) and (80) gives 


p»+! = Ti P^ 1 + l\ Q?+' (81) 

Q'lf = r 2 P'[ l + l + T 2 Q"+' (82) 

where ri )2 and Ti i2 are reflection and transmission coefficients given by 

r _ / c 2 e- 2 - cieA 

Vc 2 e 2 + ciei/ 

m _ 2e 2 ci 

— ; 

c 2 e 2 + C\C\ 

P 2 = 

\c 2 e 2 + cjei / 

rr 2e l c 2 

1 i _ ; 

C 2 e 2 + ClCl 

From (81), it is clear that a right -going wave in material 2 is a sum of a transmitted portion of a right-going 
wave in material 1 plus a reflected portion of a left-going wave in material 2. A similar argument can be 
made for the left-going wave in material 1 . In fact, the reflection coefficients ri. 2 can be shown to be identi- 
cal to the classical Fresnel reflection coefficients. The transmission coefficients also have the same form as 
the Fresnel transmission coefficients. 

Special care needs to be taken when the LBS calculates the solution at grid points near a material 
discontinuity. For example, for the x interlace at grid point (i,j) as in Figure 3, care must be exercised 
to update the solutions at grid points (i - 1 ,j) and ( i + 1 ,j). At grid point (i - 1 ,j), the term Q* +1 ■ 
in (38) becomes Q” ? -. At grid point the terms P- l and in (37) and (38) become an d Q 2:l , 
respectively. At grid point (i + 1 ,j), the term Pf-i j in (37) becomes Ff j. Rearranging equations (29) and 
(30) for grid point i we have 

(1 + 0l At) PPf - IfJi'j Id 2^) (Pfr - P ”_ ld ) - h AtQij (87) 

(1 + 02 At) Qlf - - (1 - 21*) (Q'Vl, - Qh) - b * Atp ii ( 88 ) 

where v\ — ciAt/Ax, and vo — c 2 At/Aa;. The terms a\, a 2 , b\, 6 2 refer to the a and b coefficients in (19) 
and (20) for materials 1 and 2, respectively. These equations are now easily solved for If J 1 and Q^ 1 and 
then (81) and (82) are applied to obtain F%‘^ 1 and Q 2 ,j ■ A similar analysis can be made for the boundary 
perpendicular to the y axis involving the R and S field components. 

6 Outer Boundary Condition 

The outer radiation boundary condition is used to terminate the computational lattice and permit out- 
going waves to pass unreflected through the lattice boundaries [44], The FDTD algorithm uses a spatial 
central difference operator where it uses field values from neighboring cells to update solution variables. 
Thus it cannot be used at the terminating faces of the problem domain. For example, the solution for a wave 
propagating left to right will eventually require a grid point outside the domain. To terminate the computa- 
tional lattice, an additional equation (boundary condition) is needed to solve the system and this introduces 


(83) 

(84) 

(85) 

(86) 
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information into the solution that is not required by Maxell’s equations. The PML boundary condition [45] 
has recently been introduced, which has greatly increased the accuracy of FDTD simulations. However, the 
PML comes with a moderate increase in complexity for an FDTD code due to additional variable storage 
and update equations. 

On the contrary, the LBS requires no extraneous boundary condition, and it includes the PML boundary 
condition with no extra required storage or update equations. For the present LBS implementation, like the 
Method of Characteristics [3 1 ], the interior point algorithm calculates the left -going characteristic at the left 
boundary (i.e. i — 0) and the right-going characteristic at the right boundary (i.e. i — imax ). Thus for 
the LBS, at grid point i — 0, equation (34) calculates Q(0, j) and the incoming right-going characteristic, 
P(0, j), is specified as a boundary condition. This same analysis applies at the right boundary where (33) 
calculates P(imax,j) and the incoming left-going characteristic, Q(imax,j), is specified as a boundary 
condition. Shang [20] has noted for characteristic based multidimensional and nonuniform grid problems, 
in principle, the local coordinate system can be rotated to align with the characteristics, and the compat- 
ibility equations provide an exact boundary condition. This transformation has not been implemented in 
the present work, and will likely be the subject of future studies. A simple, yet effective approximation for 
multidimensional characteristic based approaches is to set the incoming flux or characteristic variables at the 
outer boundaries to zero and let the interior point algorithm predict the outgoing variables. When the wave 
motion is aligned with a coordinate axis, this boundary condition is exact. But this approximation may not 
be necessary since the LBS automatically includes the PML boundary condition without additional storage 
or update equations. 


The linear bicharacteristic form of Maxwell’s equations for the 2D TM polarization in free space are 


dP_ dP_^ dH x 
dt + ° dx + dy 
dQ_ dQ dH x 
dt C dx ' dy 
dR dR dH y 
dt C dy dx 
d^_ dS dH u 
dt C dy dx 


0 

0 

0 

0 


(89) 

(90) 

(91) 

(92) 


In the frequency domain using complex coordinates, we have 



(93) 

(94) 

(95) 

(96) 


To show how the LBS automatically includes the PML boundary condition, we derive the appropriate up- 
date equations using the complex coordinate transformation approach proposed by Chew and Weedon [46], 
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Specifically, we use 


(97) 

(98) 

(99) 

(100) 

= 0 ( 101 ) 

- 0 (102) 

= 0 (103) 

= 0 (104) 

where B x — ( s x /s y ) H x and B y — ( s y /s x ) H y . In typical fashion with a PML FDTD implementation, we 
let <j x - a y — a, then we have that B x — H x , B y — H y and (101 )-( 1 04) become 


dP 

dt 

dP ct 

+ c— + -P + 
dx eo 

dH x 

dy 

- 0 

(105) 

dQ 

dQ a 

-c— + — Q + 
dx eo 

dH x 

= 0 

(106) 

dt 

dy 

dR 

OR a „ 

dHy 

= 0 

(107) 

dt 

t c „ ft 

dy e 0 

dx 

ds 

dt 

dS a „ 
- C — + -S- 
dy e 0 

dHy 

dx 

= 0 

(108) 


Furthermore, if we let e = eo, M = A'o and cr*//*o — <r/eo as required by the PML boundary condition, 
then the normal LBS update equations given by (2 1 )-(24 ) can easily be shown to be identical to the LBS 
PML update equations (105)-(108). This analysis shows how the LBS inherently incorporates the PML 
boundary condition within the standard update equations. The PML conductivity a is still specified using 
the conventional profiles: linear, quadratic or geometric [43]. 

7 Computational Requirements 

It is instructive to examine the computational requirements of the LBS and the FDTD method. We can 
use this analysis to determine if the LBS can provide equivalent or better accuracy than FDTD for the same 
amount of computational resources. Let us assume a 2D grid with N x N cells. The FDTD method requires 

S F - 127V 2 + 16./V + 4 (109) 


Substituting these into (93)-(96) gives 


d 

1 d 

dx 

s x dx 

d 

1 d 

dy 

Sydy 


= 1 + ■ 


jue 0 


= 1 + 


jue 0 


jojP+ — 1 
Co 

juQ+^C, 

eo 

jujR+^-1 

eo 

ju S + — . 
£o 


dP 

dB x 

dx + 

dy 

dQ 

dB x 

C dx 

dy 

dR 

d By 

° dy 

dx 

dS 

dBy 


dy dx 
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total bytes to store the field component arrays, and the LBS requires 


S L - 32 TV 2 + 24 TV 


(110) 


total bytes. Note that this storage calculation does not account for any extra terms such as arrays for boundary 
conditions, far-field transformations, etc. We can define a storage ratio §■ between the LBS and FDTD as 


S L _ 32 TV 2 + 24 TV 

~S^- ~ 127V' 2 + 167V + 4 


(111) 


If the LBS is more accurate than FDTD, we should be able to increase the cell size by a certain factor and still 
maintain the same accuracy as FDTD. Increasing the cell size decreases the total number of cells required 
in the grid. Thus, we define a grid reduction factor N r , which can be used to determine the breakeven point 
in storage and accuracy. The grid size for the LBS will be reduced in each dimension by TV, giving a new 
ratio 


32 (7V/7V r ) 2 + 24 (7V/7V,.) = 

127V 2 + 167V + 4 TV 2 *’' 


(112) 


The percentage reduction in grid storage ratio from the FDTD method is then given by 


P r = 100 '^ r - = 100 



(113) 


To determine the breakeven point, we solve P, = 0 for N r in terms of TV to yield 


7V r 


27V (4 TV + 3) 
37V 2 -I 47V + 1 


(114) 


Taking the limit of the positive root as TV -5- oo gives N r ss 1.63. Thus, the LBS must be at least 1 .63 times 
more accurate than FDTD to achieve equivalent storage for the same accuracy. Factors above 1 .63 means 
the LBS requires less storage than FDTD for the same accuracy. Figure 4 shows a plot of the breakeven 



Figure 4: Breakeven ratio versus number of grid cells, 
ratio versus the number of grid cells. 
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8 Fourier Analysis 


Various Fourier analyses of the two-dimensional LBS have already been completed [36], [37], [38]; 
therefore, only the important results and conclusions from these previous analyses will be reviewed in this 
report. Most of the information presented is summarized from [36]. The stability condition for the 2D LBS 
is v x , v y < 1/2, where v x , v y are the Courant numbers v x — cAt/Ax and v y — cAt/Ay. Although this 
stability limit is more restrictive than the standard FDTD method, it is not particularly troublesome because 
many FDTD simulations use a Courant number of 1/2 for improved accuracy. 

The complete Fourier analysis will not be outlined here for the sake of brevity. Rather, we present an 
overview of the procedure followed by a discussion numerical results. The procedure for the Fourier analysis 
is straightforward. Start with the LBS free space update equations (21)-(24) with a — b — 0 and substitute 
a solution of the form 


pn _ Pq gj (n<t>—id x —j Oy ) 


(115) 


into these expressions. After some algebra, we have the system of equations 

T n+l — Vi T n + V 2 T n ~ 1 


(116) 


which represents the three time-level LBS scheme with T" = . To complete the 

Fourier analysis, we make the substitution 'F' l+1 = T" to give 


" T ‘ 

n+ L 

’ Vi 

V 2 ' 


’ T ‘ 



h 

0 


'T 


(117) 


where h is the 4 x 4 identity matrix. The stability matrix G is then given by 


G 


Vi V 2 

h o 


(118) 


which is an 8 x 8 matrix. The stability analysis is completed by calculating the eigenvalues of the stability 
matrix G for various grid resolutions and grid propagation angles. To that end, we define 


Ox 

— 0 cos a 

(119) 

Oy 

— 0 sin a 

(120] 

0 

= 2it/N 

(121) 

4> 

= v6 

(122) 


where N is the grid resolution in cells/wavelength and a is the grid propagation angle. To simplify the 
analysis, we also set v — v x — v y . The dispersion relation can be obtained by solving the equation 

dct [e^ - G] = 0 (123) 

for <j>. In comparison, the dispersion relation for the FDTD method is 



For the one-dimensional LBS [47], it was shown the LBS had less numerical dispersion than the FDTD 
method. Extensive three-parameter studies of numerical dispersion for the 2D LBS were performed using 
the grid resolution (N), Courant number (v) and grid propagation angle (a) as parameters. These studies 
revealed that the optimum Courant number is v — 1/2 since dispersion is minimized for all propagation 
angles when compared to FDTD. 

For a Courant number v — 0.4 and propagation angle of 45°, the numerical dispersion decreases 
smoothly with increasing grid resolution as shown in Figure 5. From this figure, we see that the LBS has 



Figure 5: Phase speed error versus grid resolution N for FDTD method and LBS with v — 0.4 and a — 41? . 

approximately 1/2 the phase error as FDTD. Generally, the dispersion error for the LBS grows as v -» 0. 
When v — 1/2, numerical dispersion is zero along the coordinate axes and is maximum at 41? as shown 
in Figure 6 for a grid resolution N — 10 cells/A. When v < 1/2, dispersion for the LBS remains substan- 
tially less than for FDTD as shown in Figure 7 for N — 20. From Figures 6 and 7, it is clear that as the 
grid resolution is doubled, the numerical dispersion decreased by a factor of four; as expected for a second 
order method. Finally, as shown in Figure 8 for N — 10 cells/A, numerical dispersion decreases linearly as 
v — y 1/2; except for grid propagation angles along 45° vectors, where the LBS dispersion is very close to 
that of FDTD. For propagation along 45° vectors, LBS numerical dispersion is minimized around u = 0.3 
and then approaches the FDTD value for v = 1/2 as shown in Figure 9 for N — 20. 

To summarize, the optimal Courant number for the LBS is 1/2. This Courant number offers much lower 
dispersion for most all propagation angles except those near a 41? vector. For v < 1/2, numerical dispersion 
decreases as both grid resolution and Courant number are increased. Typically, LBS dispersion is at least 
1/2 that of FDTD, and can be much lower in many instances. 


16 





Figure 6: Phase speed error versus grid propagation angle a for FDTD method and LBS with v — 1/2 and 
N — 10. 



Figure 7: Phase speed error versus grid propagation angle a for FDTD method and LBS with v — 0.4 and 
N = 20. 
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Figure 8: Phase speed error versus Courant number v for FDTD method and LBS with N — 10 and a — CP. 
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Figure 9: Phase speed error versus Courant number v for FDTD method and LBS with N — 20 and a — 4£P . 


18 




9 Results 


To demonstrate the 2D LBS, we consider various canonical problems using the TM polarization. First, 
we inject an incoming plane wave on the outer boundaries using the LBS, and let the algorithm propagate 
the signal through the grid using a total field formulation. This is done by specifying the incoming charac- 
teristic variable (P, Q, R or S) on the appropriate outer boundary. For example, on the left x boundary, P is 
specified for all j coordinates at i — 1. We use a 71 x 71 free space grid, with a Ax — Ay = 1 cm, which 
has a time step of At — 16.67 ps and the incident wave is a Gaussian pulse with FWHM of 35 time steps 
(or 0.58 ns). We specify the incidence angle as 180% and the electric field after 160 time steps is shown in 
Figure 10. Similar results can be obtained with other incidence angles. It is clear that the LBS easily allows 



Figure 10: Propagating plane wave injected on outer grid boundaries at 18CP incidence, 
specification of incoming plane waves in its fundamental algorithm. 

Next we move on to radiation from a point source in free space. This problem demonstrates that the 
algorithm can easily treat spherical waves and it also tests the PML boundary condition. Two concurrent 
grids are used in this problem, each having a cell size of 1 mm. The first is a small test grid of size 101 x 
101 cells with an additional 10 cell PML boundary condition. This grid is centered within a large 501 x 
501 grid, and the point source is located at the center of both computational grids. The time step is 3.3 ps, 
and an electric field point source is located at the center of both grids and the total number of time steps 
is truncated at 512, to allow no reflection from the large grid outer boundaries to reach the field sampling 
points. The inner grid is terminated with PML for both FDTD and the LBS, and the large grid is terminated 
with a second-order Liao boundary condition for FDTD and a characteristic based boundary condition for 
the LBS. The electric field is sampled at the same two locations in both grids, which are located 30 cells 
in the +x direction from the point source and then ±30 cells in the y direction in the smaller grid. The 
point source is located in the smaller grid at grid point (61, 61) and the two sample points are (61, 91) and 
(61,31). Figure 1 1 shows the electric field at the upper sample point in the large grid for point source 
radiation in free space. Note the agreement is excellent, and there are no reflections from the outer boundary 
due to the Liao boundary condition. Similar results were observed at the lower sample point. Figure 12 
shows the electric field at the upper sample point (61,91) in the small test grid using the PML boundary 
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Figure 1 1 : Electric field versus time sampled at upper grid point in the large grid. 

condition. Note again the agreement is excellent. Furthermore, we computed the global error in the small 
test grid with the expression 

GE — {Elargeih j) ~ E sma u(i, j)) (125) 

ki 

using the difference between the electric fields in the large and small grids. Figure 13 shows this global error 
using the PML boundary condition for both methods and we see that the PML works very well. The error 
for the LBS is in the -80 to -100 dB range, which is excellent. Figure 14 shows the time-domain results for 
the LBS with and without the PML boundary condition. Note the reflections from the outer boundary are 
clearly visible for the no PML case. 

10 Conclusions 

This report has extended the Linear Bicharacteristic Scheme for computational electromagnetics to the 
two-dimensional case. Treatment of lossy dielectric and magnetic materials was discussed, and implemen- 
tation of the PML boundary condition was outlined. It was demonstrated that the LBS has several distinct 
advantages over conventional FDTD algorithms. First, the LBS is a second-order accurate algorithm which 
is about 2-3 times as economical. The LBS can also be made to have zero dispersion error in certain 
instances. Second, the LBS provides a more natural and flexible way to implement surface boundary condi- 
tions and outer radiation boundary conditions by using characteristics and an upwind bias technique popular 
in fluid dynamics. Third, the LBS can provide more flexibility to implement subgridding algorithms be- 
cause of the compact nature of the computational stencil. A dielectric surface boundary condition was also 
implemented and results were provided for two-dimensional free space radiation problems. Due to project 
and time limitations, validation for lossy dielectric materials and heterogeneous materials was not explored 
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Figure 12: Electric field versus time sampled at grid point (61, 91) for point source at grid point (61, 61) in 
the small grid. 



Figure 13: Global error versus time in small grid for FDTD method and LBS. 
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Figure 14: Electric field versus time sampled at upper sample point in small grid for LBS with and without 
PML boundary condition. 

in the present work. It is anticipated this will be the subject of future reports and articles. The results indi- 
cate that the LBS is a very promising alternative to a conventional FDTD algorithm for many applications. 
Higher-order extensions are available for the 2D case, but were not explored presently [36]. Extensions to 
three-dimensional problems should be straightforward. 
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